FOV Registration over Multiple Days

Progress Report Sept. 2019

In [37]:
plt.figure(figsize=[5,5])
# r, g, and b are 512x512 float arrays with values >= 0 and < 1.
from PIL import Image
rChannel = out_previous[-3]
gChanel = out_previous[-2]
bChannel = out_previous[-1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.64)*256
rgbArray[..., 1] = gChanel*(gChanel>0.64)*256
rgbArray[..., 2] = bChannel*(bChannel>0.64)*256
img = Image.fromarray(rgbArray)
plt.imshow(rgbArray)
#plt.imshow(,alpha=0.1,cmap='Reds')
Out[37]:
<matplotlib.image.AxesImage at 0x7f769a987fd0>
In [464]:
f,(ax,bx) = plt.subplots(1,2,figsize=[40,80],sharex=True,sharey=True)
A = cnmr.get_A_sparse(resultList[0])
contours = cnmr.get_contours(A,nrgthr=0.75)
ax.imshow(stack[0])
#bx.imshow(stack[0].T)
for cont in contours:
    ax.plot(cont[:,0],cont[:,1],'k',lw=1)
cm.utils.visualization.plot_contours(A,stack[0],thr_method='nrg',maxthr=0.32,nrgthr=0.5,ax=bx);
#bx.set_xlim(100,150)
In [448]:
f,(ax,bx,cx) = plt.subplots(1,3,figsize=[40,60],sharex=True,sharey=True)
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
# register each frame to the previous (already registered) one
# this is what the original StackReg ImageJ plugin uses
TrMatrix = sr.register_stack(stack[5:7], reference='first')
out_previous = sr.register_transform_stack(stack[5:7], reference='first')
################################################
RefIdx = 0
bkgIdx = 1
ARef = cnmr.get_A_sparse(resultList[5])
contours = cnmr.get_contours(ARef,nrgthr=0.85)
cont = contours[1254]
cellContour = np.append(cont,np.zeros([cont.shape[0],1]),1).reshape(-1,3)
######################################################
ax.imshow(out_previous[RefIdx],cmap='Greys_r',)
bx.imshow(out_previous[bkgIdx],cmap='Greys_r',)
for _ii,Mx in enumerate(TrMatrix[:]):
    newCoord = np.dot(cellContour,TrMatrix[0])
    for _jj in range(_ii):
        newCoord = np.dot(newCoord,TrMatrix[_jj])
    newCont = np.dot(cellContour,Mx)
    #print(newCoord)
        #plt.plot(point[0],point[1],'k.')
    #bx.plot(newCoord[:,0],newCoord[:,1],lw=1)
    if _ii==bkgIdx:
        bx.plot(newCoord[:,0],newCoord[:,1],lw=2,c='r')
        bx.plot(newCont[:,0],newCont[:,1],'k',lw=1)
    if _ii==RefIdx:
        ax.plot(cont[:,0],cont[:,1],lw=2,c='r')
        ax.plot(newCoord[:,0],newCoord[:,1],lw=2,c='k')
#ax.set_xlim(100,150)
#ax.set_ylim(250,315)
rChannel = out_previous[RefIdx]
gChanel = out_previous[bkgIdx]
#bChannel = out_previous[-1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
#rgbArray[..., 2] = bChannel*(bChannel>0.64)*256
img = Image.fromarray(rgbArray)
cx.imshow(rgbArray)
Out[448]:
<matplotlib.image.AxesImage at 0x7f76966bd898>
<Figure size 1440x1440 with 0 Axes>
In [590]:
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
f,axarr = plt.subplots(3,(len(stack)-1),figsize=[80,21],sharex=True,sharey=True)
for _ii,item in enumerate(stack[:-1]):
    out_previous = sr.register_transform_stack(stack[_ii:_ii+2], reference='first')
    ax = axarr[0,_ii]
    bx = axarr[1,_ii]
    cx = axarr[2,_ii]
    #print(len(out_previous))
    ax.imshow(out_previous[0])
    bx.imshow(out_previous[1])
    rChannel = out_previous[0]
    gChanel = out_previous[1]
    rgbArray = np.zeros((315,315,3), 'uint8')
    rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
    rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
    img = Image.fromarray(rgbArray)
    cx.imshow(rgbArray)
    ax.set_aspect(1.0)
#axarr[1,6].set_aspect(1.0)
In [588]:
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
RegStack = sr.register_transform_stack(stack, reference='previous')
f,axarr = plt.subplots(2,(len(stack)-1),figsize=[80,14],sharex=True,sharey=True)
for _ii,item in enumerate(stack[:-1]):
    out_previous = sr.register_transform_stack(stack[_ii:_ii+2], reference='first')
    ax = axarr[0,_ii]
    bx = axarr[1,_ii]
    rChannel = out_previous[0]
    gChanel = out_previous[1]
    rgbArray = np.zeros((315,315,3), 'uint8')
    rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
    rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
    img = Image.fromarray(rgbArray)
    ax.imshow(rgbArray)
    ###############################
    rChannel = RegStack[_ii]
    gChanel = RegStack[_ii+1]
    rgbArray = np.zeros((315,315,3), 'uint8')
    rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
    rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
    img = Image.fromarray(rgbArray)
    bx.imshow(rgbArray)
    ax.set_aspect(1.0)
#axarr[1,6].set_aspect(1.0)
In [592]:
baseFolder = '/media/alireza_chenani/dataWork/stress_project/315_stress_week/BL6-111/'
tiffList = sorted([item for item in locate('C*.tif',baseFolder)])
stack = np.array([io.imread(tif).T for tif in tiffList])
RegStack = sr.register_transform_stack(stack, reference='previous')
Randsin = np.random.rand(500)/2
for _ii,img in enumerate(RegStack):
    prevImg = RegStack[_ii-1]
    temp1 = prevImg*(prevImg>np.quantile(prevImg,0.9))
    temp2 = img*(img>np.quantile(img,0.9))
    overlap =  np.nansum((temp1*temp2))/np.nansum((temp1*temp1))#np.nansum((temp1*temp2)/(temp1*temp1))
    ov = []
    for sin in Randsin:
        R = np.array([
        [1-sin**2, -1*sin, 0],
        [sin,    1-sin**2, 0],
        [0,     0,         1]])
        Rimg = sr.transform(prevImg,tmat=R)
        tempR = Rimg*(Rimg>np.quantile(prevImg,0.9))
        ov.append(np.nansum((temp1*tempR))/np.nansum((temp1*temp1)))
    Thr = np.quantile(ov,0.95)
    if _ii<8:
        CLR = 'k'
    else:
        CLR = 'r'
    plt.scatter(_ii,overlap/Thr,c=CLR)
#plt.xlim(.5,14.5)
#plt.ylim(.9,8012.5)
plt.axhline(1)
Out[592]:
<matplotlib.lines.Line2D at 0x7f76913acdd8>
In [587]:
Randsin = np.random.rand(500)/2
img = stack[0]
ov= []
for sin in Randsin:
    R = np.array([
        [1-sin**2, -1*sin, 0],
        [sin,    1-sin**2, 0],
        [0,     0,         1]])
    Rimg = sr.transform(img,tmat=R)
    temp2 = Rimg*(Rimg>np.quantile(img,0.9))
    temp1 = img*(img>np.quantile(img,0.9))
    overlap = np.nansum((temp1*temp2))/np.nansum((temp1*temp1))
    ov.append(overlap)
print(np.quantile(ov,0.95))
0.3998858657661944
In [550]:
sns.distplot(ov)
plt.axvline(np.quantile(ov,0.9))
Out[550]:
<matplotlib.lines.Line2D at 0x7f769cd45630>
In [548]:
np.nansum((temp1*temp2))/np.nansum((temp1*temp1))
Out[548]:
0.7806665
In [546]:
315*315
Out[546]:
99225
In [ ]: